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Abstract 

We propose a mathematical model allowing for the alternating pulse and surge pattern of GnRH 
(Gonadotropin Releasing Hormone) secretion. The model is based on the coupling between two systems 
running on different time scales. The faster system corresponds to the average activity of GnRH neu- 
rons, while the slower one corresponds to the average activity of regulatory neurons. The analysis of 
the slow/fast dynamics exhibited within and between both systems allows to explain the different pat- 
terns (slow oscillations, fast oscillations and periodical surge) of GnRH secretion. Specifications on the 
model parameter values are derived from physiological knowledge in terms of amplitude, frequency and 
plateau length of oscillations. The behavior of the model is finally illustrated by numerical simulations 
reproducing natural ovarian cycles and either direct or indirect actions of ovarian steroids on GnRH 
secretion. 

Key words coupled oscillators, hysteresis, fast-slow dynamics, amplitude and frequency control, ovulation, 
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1 Introduction 

1.1 Endocrine background 

The reproductive function involves tightly and finely controlled processes. The reproductive axis, usually 
called the gonadotrope axis, includes the hypothalamus, within the central nervous system, the pituitary 
gland and the gonads (ovaries in females, testes in males). Specific hypothalamic neurons secrete the go- 
nadotropin releasing hormone (GnRH) in a pulsatile manner. The pulsatile GnRH secretion pattern ensues 
from the synchronization of the secretory activity of individual GnRH neurons. The release of GnRH into 
the pituitarian portal blood induces the secretion of the luteinizing hormone (LH) and follicle stimulating 
hormone (FSH) by the pituitary gland. The changes in the frequency of GnRH pulses (between 1 pulse per 
hour and 1 pulse every 6 hours in the course of an ovarian cycle) has a fundamental role in the differential 
control of the secretion of both gonadotropins: the secretion of LH is enhanced by higher frequencies, while 
that of FSH is enhanced by lower frequencies ^21- On the gonadic level, FSH and LH sustain germ cell 
production and hormone secretion. In turn, hormones secreted by the gonads modulate the secretion of 
GnRH, LH and FSH within entangled feedback loops. 

In females, the frequency of GnRH pulses is subject to the control exerted by ovarian steroids, estradiol 
and progesterone. Progesterone slows down GnRH pulsatility which leads to a slower frequency during 
the luteal phase (when progesterone levels are high) compared to the follicular phase [J]. On the contrary, 
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estradiol speeds up GnRH pulse frequency, but at the expense of a decrease in pulse amplitude so that 
the whole feedback effect on rate secretion is rather inhibitory (negative estradiol feedback). 
The GnRH secretion pattern alters dramatically once per ovarian cycle, resulting in the GnRH surge char- 
acterized by massive continuous release of GnRH in response to increasing levels of estradiol ^U] (positive 
estradiol feedback). The GnRH surge triggers the LH surge which is responsible for ovulation, leading to 
the release of fcrtilizable oocytes. 

The estradiol signal is conveyed to GnRH neurons by regulatory neurons (also designed as interneurons) 
either directly or after other neuronal relays. Transmission from interneurons calls to many different neuro- 
transmitters (see review in |51ll5)). The balance between stimulatory and inhibitory signals emanating from 
interneurons controls the behavior of the GnRH network |S] . 

Rencently, a specially interesting type of regulatory neurons has been discovered. Kisspeptin neurons act on 
GnRH neurons via the G-protein coupled receptor GPR54 They are very good candidates for relaying 
both positive and negative estradiol feedback, since they react to estradiol in opposite ways according to the 
anatomic area of the hypothalamus where they lie. 



1.2 Model objectives 

We aim at formulating a phcnomenological, data-driven model of GnRH secretion. This paper focuses on 
the coupling between the GnRH neuron network and the regulatory interneuron network. Each network is 
represented by the behavior of a single average neuron. The key point in the design of the model consists in 
entering as reasonable input 2 coupled systems (a slow one with a faster one) to generate a definite sequence 
of events in the model output: GnRH secretion. This coupling yields a 3 time-scaled model which is able 
to capture, not only the cyclic transition from a pulsatile to a surge secretion pattern of GnRH, but also 
the increase in the pulsatility frequency between the luteal and follicular phases. It also separates a specific 
dynamical state corresponding to pulsatility resumption after the surge. Besides, parametrization of the 
model is subject to physiological specifications expressed as constraints on the GnRH output and allows 
to reproduce the direct (on the GnRH network) or indirect (via the regulatory network) effects of ovarian 
steroid hormones (estradiol and progesterone) on GnRH secretion. 

In summary, we aim at reproducing a synthetic mathematical representation of GnRH secretion pat- 
tern, fitting available observations -in agreement with schematic "hand-made" representations such as that 
proposed by Herbison (top of Figure 4 in [S|). 

The paper is organized as follows. In section |21 we introduce the equations of the model. We comment the 
numerical simulations in section |2| to motivate the analysis of the bifurcations in a model-derived 2-scalcd 
system in section 0] We finally address the question of amplitude and frequency control in section |S1 



2 Model design and analysis 

2.1 Coupling oscillatory neuronal dynamics 

We consider the following four-dimensional dynamical system: 

eSx = —y + f{x) (la) 

ey = af)X + aiy + a2 + cX (lb) 

ejX = -Y + F{X) (Ic) 

Y = boX + biY + b2 (Id) 

z{t) = X{y{t)>ys} 

Equations Hla|) and (|lb|) correspond to a fast system representing an average GnRH neuron, while equations 
(|lc|l and (|ld|) correspond to the slower system representing an average regulatory neuron. The x, X variables 
represent the neuron electrical activities (action potential), while the y,Y variables relate to ionic and 
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secretory dynamics. The fast variables are assumed to have two stable stationary points separated by a 
saddle. Their bistability is accounted for by the cubic functions f{x) and F{x). The intrinsic dynamics of 
the slow variables follows a growth law of very small velocity (oi <C 1 ). In each system, the fast and slow 
variables feedback on each other. The coupling between both systems is mediated through the unilateral 
influence of the slow regulatory neurons onto the fast GnRH ones {cX term in equation ljlb|l 'l. The coupling 
term aggregates the global balance between inhibitory and stimulatory neuronal inputs onto the GnRH 
neurons. The global system exhibits 3 time scales given by e5, e and 1. Constant 7 is close to 1. 
In many cell types, thresholding calcium dynamics is known to trigger secretion. As far as GnRH secretion 
is concerned, many evidence for the inducing effect of calcium have also been gathered (see discussion in 
for details). Hence, we associate GnRH secretion to a threshold ionic activity y^i/s, and finally keep as 
representative GnRH signal z{t) = X{y{t)>y^}- 



2.2 Mechanisms underlying the pulse to surge transition 
2.2.1 Bifurcations in the fast GnRH system 

System Hla|l - (|ld| l can be analyzed in the general setting of fast-slow dynamics. The slow variable X then 
enters the fast svstem Hall bl as a parameter. Bifurcations may thus arise as the X parameter varies. The 
fast system (after changing time t into et) also exhibits fast-slow dynamics features due to the S time scale: 

Sx = -y + fix) (2a) 
y = aox + aiy + a2 + cX (2b) 
with f{x) — I'ox'^ + v\x^ -|- 1^22; (2c) 

The fast nuUcline i = is a cubic function and the stationary points are obtained as intersections of this cubic 
nuUclinc with the other nuUcline y = 0, a straight line which moves in the plane depending on the values of 
the slow variable X. Without loss of generality, we assume that lim^r^-oo /(a^) = +00. We also assume that 
the parameters (ao, a\) are chosen in such a way that there is one single stationary point whatever the value 
of X may be (see Figure ^1 . Intersection points of the cubic nuUcline with the line nuUcline are solutions of 
the equation: 

/(x) + ^. + ^^^±^ = 0. (3) 
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Figure 1: NuUclincs of the fast system 
The pink solide line corresponds to the x = Q nullelme and the cyan solid line to the y = G nullclme for X = Q (no 
coupling). The green solid line corresponds to the curve f'{x) — 0. The (xq,Xq) roots of f'{x) — — ai5~0 roughly 
correspond to the intersection points of the green line with the black horizontal line. As the value of X varies, the 
cyan nullclme sweeps the x-axis within or outside the [xq,Xq] interval, between extremal positions delimited by the 
brown and blue straight lines. 



Let xq = s(X) denote the unique solution to equation We now derive from classical arguments the 
nature of the stationary point {s{X), f{s{X))). The eigenvalues of the linearized system are solutions of the 
characteristic polynomial: 
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]A + -[ao + ai/'(a-o)] -0. 



6 ' 5' 

If the following conditions are fulfilled, two roots are complex conjugated with a negative real part. 

7/'(.T^o)+ai <0 



y + y/ i^o) > 
— ? «i 



Let Xq and Xq denote the two roots of the equation: 



-/'(x)+ai -0. 



Then, if s{X) < Xq , the stationary point is a stable focus. If X varies in such a way that xq = s{X) crosses 



the value Xq , the system (|2all - (j2b|l undergoes a Hopf bifurcation, and the stationary point becomes unstable 
as a stable limit cycle appears. This stable limit cycle disappears when xq > Xq and the stationary point 
becomes stable again. This reasoning is illustrated on Figure ^ 
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2.2.2 Limit cycle of the slow system 

The dynamics of the X variable is such that, if the straight-line nuUcline sweeps the (x,y) phase plane (from 
left to right and right to left) periodically, the Hopf bifurcations occur themselves periodically. 

The qualitative analysis of the slow dynamics is analogous to the fast one. We assume that the parameters 
{bo,bi,b2) are fixed so that there is one single stationary point (Xq, F{Xq)), which is an unstable focus 
(allowing the slow system to displays a stable limit cycle): 

-F'{Xo) + bi > 

h + hF'iXo)>0 

^ 7 
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bi 



4^<0, 



where F{x) = fio^^ + IJ'iX^ + fJ,2X. The oscillation associated to this limit cycle is of relaxation type (as 
in van der Pol system). The dynamics along the limit cycle displays slow and fast parts alternatively (see 
Figure EJ. 
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Figure 2: Slow {X, Y) hmit cycle 
The blue line corresponds to the limit cycle in the X, Y plane. The slow parts of the cycle correspond to branches of 
the cubic nullclme represented by pink diamonds (phases 2 and 4), the fast parts to the jumps from one branch of the 
cubic to the other one (phases 1 and 3). 



2.2.3 Dynamics of the coupled system 

The different phases of the limit cycle exhibited by the slow regulatory system (see Figure |2l drives the 
global behavior of the fast GnRH system (see Figure O: 

1. the first phase of the cycle, where X increases abruptly, corresponds to the ascending part of the surge; 

2. the second phase of the cycle, where X decreases slowly, corresponds to the duration of the surge; 
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3. the third phase of the cycle, where X decreases abruptly, corresponds to the decreasing part of the 
surge; 

4. the fourth and most lasting phase of the cycle, where X increases slowly, encompasses two different 
GnRH secretion patterns: 

(a) as long as X < Xq , GnRH level is almost constant, since the fast system admits a stable steady 
state and the sweeping dynamics of the straight line y = lies in a slow phase. Hence this phase 
explains the existence of a plateau after the surge. 

(b) when X > Xq_ the pulsatility of GnRH is recovered and the pulse frequency increases with X. 



3 Numerical simulations 

The numerical values of the model parameters can be constrained by physiological specifications regarding 
the features (frequency, amplitude, plateau length) of the GnRH secretory patterns [7|- The GnRH output 
should be characterized by: 

• the plateau-length to frequency ratio for the fast oscillations; 

• the pulse-amplitude to surge-amplitude ratio; 

• the surge-frequency to pulse-frequency ratio; 

• the surge-duration to whole-ovarian-cyclc duration. 

A set of parameters subject to such constraints is displayed on Tabled and the corresponding model outputs 
are illustrated on Figures 13 and ^ 
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Table 1: Numerical values of the model parameters 
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Figure 3: Outputs from the coupled systems 
pink line: X(t), cyan line: y{t) 

The main quahtative features captured by the model consist in: 

1. the cyclic transition from a pulse to a surge secretion pattern, which occurs after a short transitory 
period and seems not to be subject to initial conditions; 

2. a delay before resumption of pulsatility; 

3. the increase in pulse frequency from the luteal (post-surge) phase to the follicular (pre-surge) phase 
(see Figure El). 

If the two first properties (pulse to surge alternating and pulsatility resumption delay) were expected from 
the study derived in section |21 the third one (frequency increase) remains unexplained at this point, even if 
it is consistent with endocrinological data. Besides, the global behavior of the system is not strictly periodic, 
since both the duration of the delay and GnRH surge amplitude are not constant. The non-strict periodicity 
is a nice feature for a deterministic model in a biological context and it also reserves further analysis. We 
now explain the mechanisms underlying those properties from the bifurcation analysis of a 2-scaled system 
derived from the model. 
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Figure 5: Decrease in interpulse in- 
Figure 4: GnRH secretion pattern tervals from the luteal to the follic- 

The secretory activity of GnRH neurons occurs for threshold ionic activity ular phase (top: from time 100 to 
■^(*) = X{y(t)>ys} time 337, bottom: from time 475 

to time 710) 



4 Bifurcation analysis of a model-derived 2-scaled system 

Our approach is adapted from the " geometrical dissection" , that has been successfully applied to several 

models in Computational Neurosciences, especially those dealing with bursting oscillations. 

In classical slow/fast systems, the slow variable is "frozen" and appears as a parameter in the study of the 

bifurcations of the fast system. In a similar way, we consider the fast 3 dimensional system with 2 time 

scales: 



Sx = -y + f[x) (4a) 
y = aQX + aiy + a2 + cX (4b) 
7^: = -Y + F{X) (4c) 

where Y acts as a varying parameter. Below, we restrict our study to the case where the other parameter 
values are fixed and close to those of Tabled This fast system breaks into an independent, ID system (0cj), 
and a 2D system (|4al4b|l forced by the ID system. 

Depending on Y value, system (Pc| displays either one of the 2 possible attracting points (denoted re- 
spectively by X-(Y) and X+{Y)), or these two attractors separated by a repulsive point (denoted by 
Xa{Y)). Saddle-node bifurcations occur when Y value equals one of the local extrema of the cubic function: 
{X,Y) = ±(2/^/3, 16/(3\/3)) = ±(1.15,3.08). The whole bifurcation analysis of the 3D system IS summa- 
rized in Table 12 
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Table 2: Bifurcations of the 3D system according to Y value 
The numbers superimposed on the vertical straight lines separating the columns correspond to bifurcations of the 
following type 

(l)-(3)-(5)-(7): saddle-node, (2)- (6): saddle-node of periodics, (4)- (8): Hopf 



We now go back to the 4D system (|la|l - (|ld|) to explain the sequence of phases listed in subsection 12 .2. 31 

• In phase 1, Y decreases and system (|4a|l - (Pc|l displays the single attractive node associated to 
corresponding to the ascending part of the surge; 

• As X+{Y) decreases slowly, the solution of the 4D system remains close to the attractive node, cor- 
responding to the duration of the surge (phase 2), until this node disappears through a saddle-node 
bifurcation. Then the solution switches to the other attractive node X^{Y), corresponding to the 
decreasing part of the surge (phase 3); 

• In phase 4, as X_ {Y) increases slowly, the solution remains close to the attractive node associated to 
X_(y), corresponding to the plateau. Eventually, this attractive node disappears into an attractive 
periodic orbit initiating the pulsatile phase; 

• As phase 1 starts again, X speeds up and the pulse frequency increases. At some point the attractive 
periodic orbit disappears into a saddle node of periodics. The solution of the 4D system then jumps 
back to the single attractive node associated to X^{Y) and recovers the ascending part of the surge. 

Hence this bifurcation analysis has shown that the alternating pulse and surge pattern is based fundamentally 
on an underlying hysteresis loop. 

5 Identification of parameter targets for steroid control 

In this section we derive some tools that will be useful in controlling the amplitude and frequency of oscil- 
lations in cither the slow regulatory system or the GnRH system, hence to mimick either indirect or direct 
effects of steroid feedback. 
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5.1 Amplitude and frequency of the slow regulatory oscillations 
5.1.1 Amplitude of the slow limit cycle 

Let us first derive a parametric expression of the F{X) cubic as a function of its cxtrcma. Denote by —a 
and P the roots of F'{X) = 0. Without loss of gencrahty, wc assume a > and /3 > and choose the cubic 
function F{X) so that it passes through the origin and takes minus unity as coefficient for the higher (cubic) 
power. 

F{X) = -X^ + fiiX'^ + with ^1 > 0, ^2 > and //o = -1 

F'{X) = -3X^ + 2/iiX + ^2 

= -3{X + a){X-P) 

= -3[X^ +X{a^ (3)~a(3] 



FiX) = 

Hence a and P are the positive roots of 



X^~ -{a ^ (3)X'^ + 3apX 



3aP = fi2 

3/32 - 2fliP - fl2 



so that 



M2 



Ml + V Ml + 3^2 



and 



/3 = 



Ml + VM? + 3M2 



To compute the amplitude of the X variable, we now seek for the points of the cubic verifying either 
F{X)=F{-a) or F{X)=F{P), and F'(X)^0 (see FigureEl. 



F{X) = F(-a) 
F{X) = F{P) 



\a^{a + 3P) ^ (X^a)^ 



X 



-P^iP + 3a) {X - Pf X 



-(a + 3/3) 



-(P + 3a) 



= 



Define 



a2 = ^{a + 3/3) and P2 = ^(/3 + 3a) 



Then F{a2) = F{-a), F(-/32) = F(/3), and the global amphtude of X is thus given by a2+P2 = 2(a + /3), 
while the fast jump amplitude (related to the surge amplitude in the GnRH system) is given by a2 + a = 



P + P2 



\{a + P). In the symetric case where a = /3, those amplitudes respectively reduce to 4a and 3a 
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Figure 6: Amplitude and period of the slow limit cycle {X, Y) 
The period of the slow {X,Y) limit cycle is computed from the time taken to go along the path A(— /32, /?2)) - 
B{-a,F{-a)) - C{a2,F{a2)) - D {l3,F{f3)), with F'{a) = F'{f3) = 0, F{a2) = F{a), F{(52) = The global 

amplitude of X corresponds to a2 + 02, the surge-related amplitude to 0.2 -\- a = (3 + (32 

5.1.2 Period of the slow limit cycle 

The frequency of the regulatory oscillations drives the frequency of the GnRH surge. The period of the slow 
cycle (i.e the time taken to move along the whole cycle) in X indeed corresponds to the period of the GnRH 
surge. 

Let us denote by T the period of the {X, Y) limit cycle. The time Tx taken to move along the slow part (ie 
along the branches of the cubic function F) of the cycle can be computed as: 



Tx - 



f(-a) /-^(^ dY 



Fi-p^) boX + biY + 62 JF{a2) boX + biY + 62 
F'{X)dX [f^ F'{X)dX 



/32 



boX + biF{X) + 62 60A + biF{X) + 62 



It is worth noticing that, within the Tx duration, the time taken to climb up the descending (right) branch 
of the cubic (integration from F{a2) to F{P)) corresponds to the duration of the surge. 
The time Ty taken to move along the quasi-horizontal fast parts (corresponding to the jumps from one 
branch of the cubic function to another) of the {X, Y) limit cycle can be computed as 

, idX , ^dX 



, Fix) - Y F{X) - Y 

Alternatively, Ty can be estimated as ~ 2ejTx, according to the ratio of time constants in equations l|lc|) 
and llTd)l . 
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The whole {X, Y) limit cycle period can be finally computed from T = Tx + Ty ■ 

5.1.3 Changes in amplitude and frequency according to a prescribed ratio 

From the previous analysis, we are able to assess numerically the frequency and amplitude of the slow {X, Y) 
limit cycle for any arbitrary parameter set. We now make a further remark. Given the system: 

-fX = F{X) - Y 
y = bQX + biY + b2, 

one can one find {7, /, 60, 61, 62} to transform it into another system, 

jX = F{X)-Y (5a) 
Y = boX + biY + b2 (5b) 

oscillating with prescribed Ai frequency and A2 amplitude ratios. 
This change of variables leads to 

7 = Ai7 
bo = A160 
bi = Xibi 

O2 = 02T- 
A2 

_ \ 2 
^0 — -^2^0 

III = A2MI 

In other words, the new system reads 

7A1X = — y + Xlfinx^ + A2/^i.T^ + ii2X 
Y = X, (^boX + biY + 

The transformation does not modify the number of intersection points between the x = and y = nuUclines, 
so that the qualitative behavior of the system is preserved. 

5.2 Targeting parameters for steroid effect 

The model allows to distinguish between two possible pathways for steroid feedback on GNRH secretion, 
either directly, by acting on the parameters of the faster system, or indirectly, by acting on those of the 
slower one. The former way is dedicated to the control of the frequency and amplitude of GnRH pulses, 
while the latter is dedicated to the control of the onset time and size of the GnRH surge. 

5.2.1 Steroid-like direct effects 

Targeting adequate model parameters in the fast system in an acute way allows to alter transiently the 
amplitude and frequency of GnRH pulses. 

Figures 13 illustrates the effect of a parameter-targeting bolus (Ai = 2,A2 = 1.5), from time 525 to 575 
leading to a frequency-increased, amplitude-decreased pulsatile regime. Such an output can be compared 
with experimental results gathered in [S]"'^. These effects consist in a reduction of GnRH pulse amplitude 
and a stimulatory action on GnRH pulse frequency, as accounted by the model results. 

^ where estradiol effects on GnRH secretory characteristics are summarized in Figure 3 and GnRH portal time series are 
displaid on Figure 2 
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Figure 7: Direct, estradiol-like effect, on portal GnRH output (left) and y output (top right and zoomed on 
bottom right). The parameter values were acutely altered in a bolus way from time 525 to 575 



Figures IHl illustrates the effect of a parameter-targeting bolus (Ai = 0.5, A2 = 0.5), from time 525 to 575 
leading to a frequency-decreased, amplitude-increased pulsatile regime. Such an output can be compared 
by experimental results gathered in ^J^. These effects consist in a stimulatory action on GnRH pulse 
amplitude and a decrease in GnRH pulse frequency, as accounted again by the model results. 
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Figure 8: Direct progesterone-like effect on portal GnRH ouput (left) and y output (top right and zoomed 
on bottom right). The parameter values were acutely altered in a bolus way from time 525 to 575 



Care should be taken in the mechanistic interpretation of those simulations. Direct effect in the models 
should be rather interpreted as acute effects than direct steroid effects on GnRH neurons, whether they 
are nuclear-initiated, through steroid receptors'^ or membrane-initiated (a recently discovered signaling way 
But at the least, the simulation correspond to physiological short-term effect implying at most a few 
neuronal delays on short time scales, in contrast to the indirect, long-term effects described below. 



5.2.2 Steroid-like indirect effects 

Targeting adequate parameters in the slow system in a chronic way allow to reproduce the known effects of 
progesterone on the surge amplitude and delay for surge onset. Figures El illustrates the effect of decreasing 

^ where progesterone effects on GnRH secretory characteristics arc summarized in Figure 4 

^GnRH neurons are endowed with type /3 estradiol receptors but do not own progesterone receptors ^1] 
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the amplitude and increasing the frequency of the oscillations in the X regulatory variable. This leads to 
a decrease in the delay between two consecutive GnRH surges as well as in the surge amplitude. Such 
combined effects mimic those that have been observed in an experimental study of the long term effect of 
progesterone priming which compared the GnRH surge after exposure or not to progesterone. In the 
absence of progesterone priming, the size of the GnRH surge was decreased and its onset shortened. 
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Figure 9: Effect of progesterone- like priming on surge onset and amplitude 
Top panel: reference situation. Bottom panel situation corresponding to the absence of progesterone priming, with 
shortened surge onset and decreased surge amplitude, pink line: X(t), cyan line: y{t), blue line: GnRH output 



6 Conclusion 

We addressed in this paper the question of how the GnRH generator switches from a pulsatile towards a 
surge secretion mode. Such a question arises on the physiological scale, when considering the behavior of 
average GnRH neurons whose synchronization is taken for granted. Zeeman et al. |17| rather tackled the 
question of the GnRH-induced LH surge on the pituitary gland level, considering the GnRH self-priming on 
gonadotroph cells as a resonance phenomenon. 

On the hypothalamic level, only the variability in the frequency of GnRH pulses (rather than its control) has 
been up to now the focus of mathematical models based on nonlinear dynamics [21 El- Our modeling approach 
is comparable to these previous ones in the sense that it also considers the effect of the average activity of 
one group of neurons on the activity of another group. But the way by which this effect is introduced differs. 
They used as external inputs an impulsion train, whereas we assume that both groups can be represented 
by the same type of equations (of FitzHugh-Nagumo type) but with different time scales. Following a 3 
time-scaled approach, we have not only managed to account for the alternating pulse and surge pattern of 
GnRH secretion, but also for the frequency increase in the pulsatile regime. We have also unraveled the 
possible existence of a pause before pulsatility resumption after the surge, which could be investigated from 
an experimental viewpoint. Hence the capacity of our model to display complex features interprctablc against 
experimental evidence suggests that such a modeling approach may be a useful complement to experimental 
studies of neuro-endocrine systems. 
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